/*
 * AnalyzeCharges.h
 *
 *  Created on: May 26, 2011
 *      Author: marchi
 */

#ifndef ANALYZECHARGES_H_
#define ANALYZECHARGES_H_

#include "gmx_charges.h"
#include <iostream>
#include "Metric.h"

static int AnalyzeCharges(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
    /* Here, you can do whatever analysis your program requires for a frame. */

	static int nframes=0;

	int ecenter=0;
    int ePBC;

	int unitcell_enum=euCompact;

    gmx_rmpbc_t  gpbc=NULL;

    //    int ru=gmx_ana_get_topconf(trj,NULL, NULL, &ePBC);
    ePBC=fr->ePBC;

	set_pbc(pbc,ePBC,fr->box);

	gpbc = gmx_rmpbc_init(&top->idef,ePBC,fr->natoms,fr->box);

    gmx_rmpbc_trxfr(gpbc,fr);

	put_molecule_com_in_box(unitcell_enum,ecenter,
				&top->mols,fr->natoms,top->atoms.atom,ePBC,fr->box,fr->x);

	reset_x(natoms, NULL, natoms, NULL,fr->x,w_lst);

    std::cout << std::endl <<  fr->x[10][0] << std::endl;
    if(!nframes){
    	copy_rvecn(fr->x,xrefs,0,natoms);
    }
    std::cout << fr->x[10][0] << std::endl;
    matrix R;
    calc_fit_R(DIM,natoms,w_lst,xrefs,fr->x,R);

    Metric Met;
    Met(fr->box);
    for(int i=0;i<ngrps;i++)
    	chg[i].setCoord(Met,fr->x,cidx[i]);

	std::cout << chg[Sol][2][0] << " " << chg[Sol][2][2] << " " << chg[Sol][2][2] << " " << std::endl;
    exit(1);

	Rho * ro=new Rho [ngrps];
	Phi * fi=new Phi [ngrps];
	unsigned int nx0=static_cast<unsigned int>(ngrid[0]);
	unsigned int ny0=static_cast<unsigned int>(ngrid[1]);
	unsigned int nz0=static_cast<unsigned int>(ngrid[2]);
	ro[Slt].Allocate(nx0,ny0,nz0);
	ro[Sol].Allocate(nx0,ny0,nz0);
	ro[Tot].Allocate(nx0,ny0,nz0);

	fi[Slt].Allocate(nx0,ny0,nz0);
	fi[Sol].Allocate(nx0,ny0,nz0);
	fi[Tot].Allocate(nx0,ny0,nz0);

	ro[Slt].Density(chg[Slt]);;
	ro[Sol].Density(chg[Sol]);;
	ro[Tot].Density(chg[Tot]);;
	std::cout << " rho slt " <<ro[Slt].Sum() << std::endl;
	std::cout << " rho Sol " <<ro[Sol].Sum() << std::endl;
	std::cout << " rho Tot " <<ro[Tot].Sum() << std::endl;

	fi[Slt](ro[Slt]);
	fi[Tot](ro[Tot]);
	fi[Sol]=fi[Tot]-fi[Slt];

	exit(1);

	fprintf(myfp," Frame No. = %d \n",nframes);

	if (gpbc != NULL)
    {
        gmx_rmpbc_done(gpbc);
    }


    /* We need to return 0 to tell that everything went OK */
	nframes++;
    return 0;
}


#endif /* ANALYZECHARGES_H_ */
